Prenet algorithm for simulations of fluctuating continuous elastic filaments 
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We present a new algorithm for generating the equilibrium configurations of fluctuating continu- 
ous elastic filaments, based on a combination of statistical mechanics and differential geometry. We 
use this to calculate the distribution function of the end-to-end distance of filaments with nonva- 
nishing spontaneous curvature and show that for small twist and large bending rigidities, there is 
an intermediate temperature range in which the filament becomes nearly completely stretched. We 
show that volume interactions can be incorporated into our algorithm, demonstrating this through 
the calculation of the effect of excluded volume on the end-to-end distance of the filament. 
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Theories and computer simulations of polymers are based on the notion that a macromolecule can be modeled 
as a collection of points with positions {rj} that represent either "real" chemically bonded monomers interacting 
through semi-empirical potentials (e.g., in Molecular Dynamics simulations [|] ) , or statistically independent segments 
connected by elastic springs or subject to constraints (in Monte Carlo simulations [Q). In the latter case, the elastic 
energy is usually assumed to be of entropic origin|^ E^i = {K/2)J2i{''^i — Tj-i)^, where K = ksT/aL is the spring 
constant, ks the Boltzmann constant, T the temperature, L the polymer length and a the monomer size. In polymer 
physics one often employs the continuum version of this theory, the so-called Gaussian chain (GC) model, in which 

the monomer label is replaced by a continuous contour parameter s, — *■ r(s) and Eqc — {K' /2) ds [dr/ds)"^ 
{K' is a constant). Alternatively, one can use the worm-like chain (WLC) model in which the energy penalty for 

stretching of elastic springs is replaced by bending elasticity, Ewlc = {kBTlp/2) jj" ds (d^r/ds^) where Ip is the 
bending persistence length. At first sight it appears that in taking the continuum limit we pass from a description of 
a polymer as a collection of points to one in which it is described as a line. However, any geometrical line in 3d space 
must obey the Pythagorean theorem, \dY{s)lds\ = 1, a condition that can not be imposed in the framework of the GC 
model (it would make Eqc a conformation- independent constant). This is consistent with the observation that the 
statistical properties of the GC model are identical to a those of a random walk and therefore the conformation of a 
polymer in this model is a fractal, with fractal dimension 2 (recall that the fractal dimension of a line is the same as 
its geometric dimension, 1). Although the above constraint is consistent with EwLC^ the resulting nonlinear theory 
appears to be intractable. Nevertheless, it was shown that the statistical mechanics of this model can be worked 
out using the analogy between the WLC and a quantum top, and the results were successfully applied to model the 
stretching of DNA molecules Q. 

Recently, we analyzed the statistical mechanics of the generalized WLC model that describes the linear elasticity 
of ribbons, with elastic energy 
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where the coefficients bi and 62 are associated with the bending rigidities with respect to the two principal axes of 
inertia of the cross section (they differ if the cross section is not circular), and 63 represents twist rigidity. In this 
paper we will treat {6^} as given material parameters of the filament. The functions {wfc(s)} and {uJko{s)} are related 
to the generalized curvatures and torsions in the deformed and the stress-free states of the filament, respectively. 
These parameters completely determine the three dimensional conformation of the centerline and the twist of the 
cross-section about this centerline, through the generalized Frenet equations 

^ ~l^£i3kl^3{s)tk{s). (2) 

Here ta is the unit tangent to the centerline and the unit vectors ti and t2 are oriented along the principal axes of the 
cross section {sijk is the antisymmetric tensor). Note that since these equations describe a pure rotation of the triad 
{ti(s)} as one moves along the contour of the filament, the constraint jtaj = \dY{s)/ds\ = 1 is automatically satisfied 
in this intrinsic coordinate description. 

Since the energy is a quadratic form in the deviations 5u)k = — w^-o of the curvature and torsion parameters from 
their values in the stress free state (Eq. (|^) is valid as long as the characteristic length scale of the deformation is 
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larger than the diameter of the filament |^ ) , the distribution of Sujk is Gaussian, with zero mean and second moments 
given by|| 

{5iu,{s)SLJ,{s')) = ^6,j6{s - s'). (3) 

Using the above expression we showed that all the two-point correlation functions {ti{s) ■tj{s')) can be calculated 
by solving a simple linear differential equation with s-dependent coefficients, for arbitrary parameters of the stress 
free state {ujko} and rigidity parameters {bk}- However, since the distribution functions of the various fluctuating 
quantities (e.g., the end-to-end distance R) are non-Gaussian, knowledge of the second moments does not determine 
the distributions and therefore, the complete determination of the statistical properties of fluctuating ribbons requires 
more powerful analytical or numerical methods. In this paper we present an efficient algorithm for the simulation of 
fluctuating elastic ribbons and use it to study the effects of spontaneous curvature and twist rigidity on the spatial 
conformations of fluctuating ribbons. To the best of our knowledge this is the first direct simulation of continuous 
lines (other simulations of the WLC model and its extensions represent the filament as a collection of points 0, 

Examination of Eq. (^ shows that the decoupling of the "noises" {Suk} in the intrinsic coordinate representation 
permits an efficient numerical generation of independent samples drawn from the exact canonical distribution. The 
Gaussian distribution of the Sujk means that each Sujk{s) can be directly generated as a string of independent random 
numbers drawn from a distribution symmetric about the origin with width y^kgT/bkds, where ds is the discretization 
step length. Note that the discretization of the continuous filament (choice of ds) is a computational tool only, and 
is very different from the case of a chain consisting of discrete mass points. We always choose ds sufficiently small so 
that the results are insensitive to its precise value. 

The remaining task is to construct the curve using the Frenet equations with ujk{s) = ujko{s) + Sujk{s)- The Frenet 
equations are best integrated by stepping the basic triad {t^} forward in s through a suitable small rotation. In this 
way, the orthonormality of the triad is guaranteed to be preserved up to machine accuracy. To construct this rotation 
matrix, we begin with Eq. (||) and, defining the three vectors — (if , tf > '•Dj and so forth, we can write this equation 
as 

^ = Av^ (4) 
ds 

where A is an antisymmetric matrix with elements Aij — SijkLOk- We now discretize Eq. (^ as 

v\s ^ ds) ^ Ov\s) (5) 

where the orthogonal matrix O is 

o.(:.|.)(,-|.)-' ,e, 

In the following we present the results of simulations of distribution of the end-to-end distance for a filament of 
contour length L = 1 and study its dependence on the stress free configuration {ajfeo} , rigidity parameters {6fc} and 
temperature. 

We first consider a straight filament, with wio = lo2q = = 0. In the limit T ^ 0, the distribution P{R) of the 
end-to-end distance, R, approaches a delta function peaked at i? = 1. With increasing T, the peak of the distribution 
shifts to smaller values of R, consistent with the fact that the effective persistence length scales as Ip ^ h/ksT [h is a 
characteristic bending rigidity parameter). The behavior of the width of the distribution is more interesting. Initially 
(in the range 3> 1) the distribution broadens with T and then narrows down again as the Gaussian chain limit 
(R^) = IpL oc 1/T is approached for Ip <C 1. This behavior is universal and takes place for arbitrary values of the 
rigidity parameters; furthermore, the form of the distribution depends only on the bending rigidities bi, 62 and is 
unaffected by the twist rigidity 63. 

Now consider a ring, broken at a point so that the ends are free to fluctuate. Here uiq — 27r, LU20 = 1^30 = 0. In 
Fig. la we present the distribution function P{R) for the case 61 = 62 = &3 = 1- As dictated by the geometry of 
the stress free state, P{R) approaches a delta function peaked at i? = in the limit T — > 0. At higher temperatures, 
fluctuations increase R, the distribution broadens and its peak moves out to higher values of R of the order of the 
diameter of the ring (I/tt). At yet higher T the decrease of the persistence length with increasing temperature takes 
over, the distribution narrows and its peak moves to smaller values of R. 

In Fig. lb we present the case of the broken ring now at small twist rigidity (63 ^ 62). While the low temperature 
behavior is similar to that in Fig la, nearly full stretching accompanied by a dramatic narrowing of the distribution 
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FIG. 1: Plot of the distribution function P{R) vs. the end-to-end to end distance R for an open ring, a) 61 = 62 = fes = 1, 
and temperatures T = 0.1 (circles), 1 (squares), 10 (crosses) and 100 (triangles), b) 61 = 62 = 1, &3 ~ lO"**, and temperatures 
T = 10~^ (triangles), 10~^ (crosses), 0.1 (pluses) and 10 (squares). Snapshots of typical configurations corresponding to each 
of these temperatures are shown as inserts. 
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FIG. 2: Simulation results for the mean square end-to-end distance {R^) vs. temperature T for an open ring, with bi — 62 — 1 
and 63 = 1 (circles), 63 = lO"'^ (squares) and 63 = 10~^ (triangles). The results of the analytical calculations are shown as 
solid lines. 



of the end-to-end distance is observed at intermediate temperatures for which ai/L — a2/L > 1, a^/L -C 1 (we 
define the bare persistence lengths, {a^} = {bi/ksT}). This striking observation is supported by analytical results 



where the elements of the matrix A 
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for the mean square end-to-end distance^ (R^) — 2 ds ds' e * ) 

are A^^ — {kBT/2) J^j ~ ^ik + J2i ^iki^io- An exact expression for (R^) in terms of the parameters {bk} 
was derived by diagonalizing the matrix A and is compared against our simulation results in Fig. 2. We would like to 
emphasize that the nearly complete stretching of the filament is a finite L effect and that the renormalized persistence 
length defined by Ip = lim^^oo (R^) /L, does not diverge in the limit of vanishing twist rigidity. 

The observation of nearly complete stretching of the filament is at first sight counterintuitive, given that at these 
temperatures, its bending persistence length is larger than its contour length, and so thermal fluctuations can not 
significantly change the curvature from its finite spontaneous value (wio = 27r), let alone cause it to vanish. The 
resolution of this paradox is that our intuition that a circle can not be continuously deformed into a line, at fixed 
curvature, is wrong. Any 3d curve can be described by its curvature wi and torsion lut, (if we replace the filament 
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FIG. 3: Simulation results for the mean square end-to-end distance {R^) vs. the contour length L for a linear filament, with 
(circles) and without (squares) excluded volume. The solid lines give the least square fit to the data. The parameters are 
6i = b2 = 63 = 0.02, and T = 1. 

by a geometrical line, we should set[|| ^2 = 0). Keeping the curvature fixed and increasing the torsion (assumed to 
be constant along the filament), one goes from a circle of radius wj"^ to an increasingly stretched helix, with radius 
uji/ (ujI -I- w|) — > and end-to-end separation approaching that of a straight line. This indeed takes place in our 
case since Eq. (^) gives (i^l) = ksT / {h'ids) so that thermal fluctuations of the torsion ^3 are large in the limit of 
small twist rigidity. Since, under conditions of nearly complete stretching, the distribution of the end-to-end distance 
is extremely narrow (see Fig. lb), mean field considerations apply and we can describe the filament as an object 

with length L = 1, curvature coi = 2tt and torsion ((^3) = ±y/kBT/{b3ds). For high temperature and small twist 

rigidity, | | wio, this object is a nearly straight helix that oscillates between left-hand and right-hand forms 

as one moves along its axis (in this limit the change of sign of torsion changes the sense of helical rotation but does not 
affect the direction of the axis about which the helix rotates), with typical radius of helical turn 2Tihj,ds/ {ksT) that 
is much smaller than its pitch l-K^h^dsj (ksT). Since in this limit the pitch approaches the period of the helix, we 
conclude that the helical filament is nearly completely stretched. At yet higher temperatures, the bending persistence 
length becomes shorter than the length of the filament, the helical structure "melts", and one recovers the Gaussian 
chain behavior characteristic of fiexible polymers. 

In order to demonstrate that volume interactions between parts of the filament can be readily incorporated into 
our algorithm, we calculated the Flory exponent ly defined by the relation, (R^) ^ L^'^ . We inserted 50 equally 
spaced interaction sites per unit length into the filament (50L sites in total) and allowed only configurations in which 
the spatial distance between any two of these sites exceeded 0.01. For each value of L we generated 1000 allowed 
configurations and calculated the average value of R^. We first checked that when all configurations were allowed 
(no excluded volume), we obtain the Gaussian random walk result, vq « 0.51. In the presence of excluded volume, 
the best fit to the simulation (Fig. 3), yields vsaw ~ 0.59, in excellent agreement with the Flory exponent of 3/5 for 
self- avoiding polymers in 3d[^. 

In closing, we would like to comment on possible ramifications of this work. The methods presented in this work 
are ideally suited to investigate polymers with spontaneous curvature such as double stranded DNA[Q, and synthetic 
molecules with bent cores |ll|. As far as equilibrium properties are concerned, our algorithm is computationally far 
superior to standard Monte-Carlo methods in that each realization of the filament is completely independent. Because 
sequence dependence can be readily incorporated into the present theory through the dependence of the spontaneous 
curvature and torsion parameters {uJko{s)} on the contour position s, and because the method is easily extendible to 
incorporate excluded volume, electrostatic and other self-interaction effects, the Frenet algorithm has the potential of 
becoming an important new tool for computer simulations of equilibrium properties of complex biopolymers such as 
DNA, RNA and proteins. 
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